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ABSTRACT 


Natural variability in mean temperatures and total precipitation 
has been related to fluctuations in the levels of lakes contained 
within the Cooking Lake Moraine area. 

The monthly lake levels have been estimated by developing a 
water-balance model from historical meteorological data recorded at 
Edmonton, physical characteristics of the basin and lake level records 
of Cooking Lake. The technique devised by Thornthwaite, et al, has been 
used to obtain the best estimate of evapotranspiration and surface 
runoff, the major components of the model. 

This thesis shows that lake level fluctuations can be related 
to natural variations in meteorological variables, by generating a 
long record of meteorological data using Monte Carlo Methods, and by 


obtaining simulated lake levels using the water-balance model developed. 
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CHAPTER I 


INTRODUCTION 


The objective of this study is to develop a mathematical model 
which will describe the stochastic nature of the water balance in 
the Cooking Lake Moraine area, located approximately eighteen miles 
south-east of Edmonton. Figure t is a map of the area showing the 
main lakes of the system and the approximate flow path of streams 
connecting the lakes. The water balance model will be developed from a 
consideration of the available data and of the variables which effects 
directly or indirectly the input or output of water contained within 
the lake system. 

An experimental statistical technique will then be used to generate 
a series of psuedo random numbers with a given distribution function. 
The given distribution function will either be the same as or within 
acceptable limits as that of the observed variable. The generated 
values will be used as the input variables of the water balance model 
and an investigation made of the expected fluctuations in the levels 
of the lakes contained within the boundaries of the Cooking Lake Moraine. 
The results of this study are applicable in a broad sense to all lakes 


in the same climatological region. 


Background Information on the Area 


In 1970 a petition, signed by local citizens of the Cooking and 
Hastings Lake area, called for Provincial Government action to reclaim 
the Beaverhill Watershed. As a direct result of this petition the 


Province organised and financed a multi-disciplined study of the entire 
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Cooking Lake Moraine to determine its potential as a recreational 

area for the City of Edmonton and as a standard for other conservation 
areas within Alberta. [Stanley Associates, 1974. EPEC Consulting, 
1971, 1976. ] 

The prime location of the area, Figure 1, makes it an invaluable 
source of water based recreational activities. This potential has 
been recognised since the first settlers moved intothe area and the first 
recreation establishment was the Goney Island Sporting Club formed in 
1894. This club was basically a rather exclusive organisation available 
to only the more affluent members of the community. It was not until 
1909, when the first passenger train of the Grand Trunk Pacific Co. 
ran passed the north shore of Cooking Lake, did the area become avail- 
able to the working classes of Edmonton. 

The area is generally known as the Cooking~-Lake Moraine but it also 
has several alternative names, principally Beaverhill Watershed area, 
and its old Indian name Amisk Wuchee (Beaver Hills). In whatever case 
it still comprises an integral series of lakes believed to be dependent 
to various extents on each other. 

The topography of the area is characteristic of hummocky disinte- 
gration moraine. This glacial moraine was formed during and after the 
last glacial period,of; about 9000 years)ago..,The level of ,the,area is 
generally above that of the surrounding land making a watershed area of 
approximately 1500 square miles, draining from the south into the North 
Saskatchewan River via Beaverhill Lake. 

Early descriptions of the lake system and its surrounding area 
mentions the abundance of pelican, cormorants and blue heron nests, 


which implies a large fish population, thousands of deer, elk and moose 
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and dense spruce forests surrounding clear lakes connected by streams 
Suitable for large canoes (Nyland, 1969). The past eighty years has 
seen this description change to one where the wildtife has been decimated, 
the forests cleared and the lakes, in certain cases, are mere sloughs. 
The Beaver Hills were invaded around 1780 by the Blackfeet during 
their migration from Eastern Canada. They, in turn, had to yield to 
the Cree. Around 1890 the influences of the first white settlers began 
to be felt. Almost immediately upon their arrival to the area the settlers 
began to clear the forests, with fires, so that cultivation of the land 
could proceed. Some of these fires burned out of control due to dry 
ground conditions, high winds or the over zealousness of the fire setters. 
In several cases large areas of prime forest were destroyed, the most 
extensive occurring in 1895. This fire destroyed an area of forest 
bounded by Edmonton in the west, Beaverhill Lake in the east, Fort 
Saskatchewan in the north and Cooking Lake in the south. 
A particularly important event, with respect to the water balance 
of the area, was the opening of a canal from Miquelon Lake to a seseredit 
serving the City of Camrose in 1927. Up until the opening, of the canal, 
the lakes were in good condition but the opening year and the first years 
of operation coincided with relatively dry years and a drastic drop in 
the level of Miquelon Lake occurred. Three shallower lakes formed 
from Miquelon Lake after this drop. The lakes downstream of Miquelon 
Lake also appear to have been affected. It has been reported (Nyland, 
1969) that in the dry years following the opening of the canal the 
water levels in Oliver, Ministik and Cooking Lakes were lower than 


anyone could remember. 


Although the canal is no longer in use and has been blocked off, 
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Miquelon Lake has been unable to recover naturally to its former regime. 
It would therefore appear that 1927 marks an obvious turning point 

in the water balance of Miquelon Lake and possibly the entire lake 
system. 

Activities such as land clearing by fire setting, water diversions 
and changes in land use have all played their respective roles in the 
decline of the lake system. It is the purpose, however, of this treatise 
to investigate the important role nature plays in the water balance of 
the system, in particular, the stochastic variations in precipitation, 
evaporation and evapotranspiration, and to determine if natural variabil- 
ity in weather could have caused the large fluctuations in lake levels 


that have been observed. 


Surficial and Bedrock Geology 

The topography associated with hummocky disintegration moraine 
characterises the Cooking Lake Moraine area. Hummocky moraine, by 
definition, has a local relief of more than fifteen feet and is generally 
thick (40 to 150 feet). [Bayrock and Hughes, 1962. Bayrock, 1972]. 
The composition of the moraine is primarily of till made up of almost 
equal proportions of sand, silt and clay containing pebbles and boulders. 
Some lenses of sand, gravel and local bedrock are also to be found in 
the surficial deposits. Recent lacustrine deposits have been mapped 
around Cooking, Joseph and Miquelon Lakes and it is believed that the 
other lakes in the system will show the same general deposits. [Bayrock, 
1972]. The lacustrine deposits are predominately composed of silt and 


clay with local marl deposits, with some clean sand found in places. 
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Several geological cross-sections and surficial features are defined 
in Figure 2 and the cross-sections shown in Figure 3. [Stanley & 
Associates, 1974]. 

The surficial deposits described above overlie . the Upper Cretaceous, 
Horseshoe Canyon formation of the Edmonton group. The Horseshoe Canyon 
formation is underlayin by the shaley beds of the Bearpaw formation 
which is in turn underlayin by the Belly or Judith River formation. The 
Belly formation has a lithology similar to that of the Edmonton group. 
[Carlson, 1966]. 

The main features of the bedrock topography are a general coincidence 
of preglacial and present day water divides, and major buried valleys 
following the general trend of Hastings and Katchemut Creek. [Carlson, 
1966] There have been two branches of the Vegreville bedrock channel 
identified in the area. Ministik lake lies over one of these branches 
and Cooking Lake over the other [Farvolden, 1963. Carlson, 1966]. Sand 
and gravel in places, exceeding thirty feet in thickness, have been 


reported within the buried valley beneath Cooking Lake. 
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FIGURE 3: Geological Cross-sections (after Stanley and Associates) 
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CHAPTER 2 
WATER BALANCE MODEL 

Introduction 

The Cooking Lake Moraine contains a system of interconnecting 
lakes which can, depending upon the lake levels, contribute to the 
next downstream neighbour. A primary flow series can be seen from 
Miquelon Lake through Larry Lake to Oliver Lake. Oliver Lake to 
Ministik Lake to Cooking Lake and finally to Hastings Lake will complete 
the primary series (Figure 4). The flow continues from Hastings Lake 
to Beaverhill Lake which discharges to the North Saskatchewan River. 
Primary Lakes are defined as those which are on the main flow systen, 
receive a contribution from another lake and contribute to other lakes. 
Secondary Lakes are those which only contribute to the water balance 
of other lakes. The secondary system consists of Joseph Lake, Half—Moon 
and Antler Lakes, and Wanison Lake. Figure 4 shows a cross-sectional 
profile of the Cooking Lake Moraine with a few reference elevation lines, 
and Figure 5 is a schematical flow chart illustrating the idea of primary 
and secondary flow systems. 

The main components of any regional water balance model are shown 
in Figure 6. The interdependence of these components can be seen by 
the inspection of the following equation (Eqn. 1) which mathematically 
defined the water balance relationship with respect to lake level 


fluctuations. 


AS = PPT + IF +GW+R- OF-E (1) 
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FIGURE 5: 


Flow Chart of Lake Dependence 
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FIGURE 6 


REPRESENTATION OF WATER-BALANCE EQUATION COMPONENTS 
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E : Evaporation from lake surfaces 


R : Runoff from land after storage capacity of ground has 
been filled 

GW: Ground water input into the lake watershed 

AS : Change in lake level or water stored in lake 

OF : Outflow from lake to the next watershed 

PPT:* Precipitation 

IF : Inflow into the watershed 


ASC: Change in ground storage 
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where: 
AS = change in lake storage 
PPT” = "precipitation on lake and land surfaces 
IF = inflow from streams and other lakes 
R = runoff from the lake's watershed 
OF = outflow from the lake 
GW = groundwater contribution 
E = evaporation from lake surface 


A convenient time interval is usually selected and is dependent 
upon the interval at which the independent variables of equation (1) 
are measured. All of the above components will vary with the time 
of year and with the physical and climatical conditions prevailing at 
any particular time. Each component in turm is considered to deter- 


mine its significance with respect to the lake level fluctuations. 


Precipitation 

As in many hydrological studies the lack of suitable measurements 
will hinder the development of a model and will necessitate the usage 
of basic simplifying assumptions, which may or may not be valid. Precipi- 
tation records for the immediate Cooking Lake area are unavailable for 
the length of period required. It is therefore necessary to adopt 
records of neighbouring stations. Edmonton and Edmonton International 
Airport are two such stations situated approximately twenty miles 
north-west and south-west respectively, from Cooking Lake. The record 


periods are, 1880 to present for Edmonton and 1960 to present for Edmonton 


International Airport. 
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Numerous investigators of hydrological phenomena have noted that, under 


particular conditions, naturally occurring variables such as precipi- 
tation and temperature conform to a normal frequency distribution. 
[Markovic, 1965. Yevjevich, 1972]. The normal distribution was probably 
first described by Abraham De Moivre in 1733 and used later by Gauss 

to describe the theory of accidental errors of measurements involved 

in the calculation of orbits of heavenly bodies. For the normal distri- 
bution only two parameters, the variance and the expected value, are 
needed to completely describe the distribution. 

The remainder of this section will be devoted to the determination 
of the frequency distribution associated with total monthly precipitation. 
When this distribution is known a sample of some length (say 500 years) 
will be generated by the experimental statistical technique (Monte Carlo 
Method). 

The normal distribution has been observed to be followed by many 
hydrological variables whenever these variables are taken in their 
natural or in some transformed state. It is therefore obvious to start 
an investigation of this type by testing the normalcy of the available 
time series. The theoretical fitting of the normal distribution to an 
empirical distribution often causes difficulties. One of the main 
problems is associated with the property of the theoretical normal 
variate being able to assume any value between plus and minus infinity, 
whereas the many hydrological variables measured, assumed to be normally 
distributed, are positive valued with a lower boundary, usually zero. 
Moreover if the properties of the normal distribution are recalled, the 
probability of the variate being negative or reaching the lower boundary 


is very small if the mean is greater than three standard deviations. 
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In many cases this minute probability can be neglected and the bounded, 
positive valued variable can be assummed to be normal. 

Several precedures are currently in mode to show the applicability 
of the normal distribution to an observed time series. The simplest 
is a frequency plot of the magnitude of the variable against a plotting 
position, usually equal to the rank of that value divided by N+ 1, 
where N equals the total sample length. The graphical method for the 
determination of the parameters of a distribution is’limited in its 
accuracy as several investigators will produce different results, but 
it does show how plausible the normal fit will be. The plots of total 
annual precipitation, natural log and cube root of total annual precipi- 
tation, mean and natural log of mean annual temperature are given in 
appendix 3 as an example of the type of plots produced. The best pro- 
cedure, to standardise results, is to use a mathematical approach to 
obtain the sample's distribution parameters. 

There are basically four methods of estimating the population 
parameters of a distribution function from a sample of that population. 
Generally these are known as, the least squares method, the method of 
moments, the maximum likelihood method and the graphical method described 
earlier. The first three effectively give the mean and variance of the 
normal function to the same accuracy. 

The population mean is estimated by the sample mean, X, and the 


population variance by the unbiased sample variance, 5 . For the case 


of a discrete variable, the mean is given by 


N 
) x (2) 
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and the variance by, 


K N sae 
unbiased, s? aye ) Coes) =) (3) 
Nee j 
j=l 
where, 
N = number of observations 
= = jth observation of variable x. 


and other terms defined above 


In terms of moments; equation (2) defines the first moment about 
the origin and equation (3) defines the second central moment about the 
mean with a correction for bias. Using equations (2) and (3) the popu- 
lation estimates of expected value and unbiased variance for total 
monthly precipitation were obtained. These estimates are given in 
table 1. Also given in table 1 are the estimates of the population's 
skewness and kurtosis coefficients. The skewness coefficient can be 
derived directly from third and second central moments and the coefficient 
of kurtosis from the fourth and the second central moments. Mathemati- 
cally the definition of skewness is given by equation (4) and kurtosis 
by equation (5). The Excess Coefficient is another useful distribution 


parameter and is defined by equation (6) 
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E = a =~ 3 (6) 
where 
on = unbiased estimate of the population skewness coefficient 
om = unbiased estimate of the population coefficient of kurtosis 
9 = unbiased estimate of the population excess coefficient 
and 
M, = ith central moment 
ae: =, 
= noe (x, ~ §x) 
j=l 


The coefficient of skewness is a measure of the asymmetry indepen- 
dent of the dimension of the variable and the coefficient of kurtosis 
measures peakedness or flatness. These two estimates of the population 
parameters provide useful tests to the validity of assumming normalcy 
of a time series. This is due to the normal distribution having a 
coefficient of skewness equal to zero and the coefficient of kurtosis 
equal to three or the excess coefficient equal to zero. 

The following two identities have been used to test for Normalcy. 
(Yevjevich, 1972). 

(1) Unbiased skewness coefficient, G. satisfies either 

(a) -0.10 < G, < 0.10 (7) 
or more confining 

(b)  -0.05 pao < 0.05 (8) 
and, (2) Unbiased Excess Coefficient 


=O. S ob Uae (9) 
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TABLE 1 
Unbiased population estimates of mean, variance,coefficient of 
skewness and coefficient of kurtosis for total monthly precipitation 


oe k 
measured in inches at Edmonton (Industrial Airport). 


Edmonton N = 91 


Month Mean S C. Cc. 
JAN 0.93 0.54 0.70 3.44 
FEB UAT Ae 0.23 eer 4.34 
MAR Ooi, O55 153 6.76 
APR 0593 0.67 1.41 0 
MAY 1.76 nee eS) Oy, JAGR sis) 
JUNE 3.14 1.56 Org 329 
JULY Si56S! (eee) eo 0) oon 
AUG Do JU) 1.42 0.54 2.69 J 
SEPT Hippos 0.93 1.00 3.299 
OCT O74 0.56 0.67 yee | 
NOV 0.74 0.56 iy S17 Pao 
DEC 0.83 Ono7 1.40 70 
ANNUAL A I Be es 3.66 WA ¥ 4.06 


* Monthly Records - Atmospheric Environment Service 
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For the ninety-one years of records recorded at Edmonton, estimates 
of both the skewness and kurtosis coefficients have been obtained and 
the identities given by equations 8 and 9 applied. The conclusion 
drawn from these two tests is that the variable, total monthly precipi- 
tation, does not have an empirical distribution function that is close 
to the normal distribution function. Fortunately there are two alter- 
natives still available. The first alternative is to try and fit some 
other theoretical distribution to the observations and the second is 
to transform the variable so as to obtain a 'new' variable which may 
be normally distributed. 

The second alternative is the more obvious and convenient path to 
follow at this stage. A transformation will still allow the use of 
the simple properties of the normal distribution and prevents the use of 
a theoretical distribution function which may give a close fit but would 
be difficult *tosjustify: 

Transformations to observed data take many forms, the most common 
being the natural logarithm. This transformation proved to be not 
normally distributed, that is, not log-normal, when the procedure and 
tests described earlier in this text were applied. The next most common 


transformation usually takes the form described by equation (10). 


it 6 ae ee (10) 


where: 
X = Random variable 
T(X) = transformed variable 
I = transforming agent. 
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For the time series, given by the records of total monthly precipi- 


tation, a series of I values were applied and the distribution parameters 


“a “nw “an 


of, expected value, variance, Co. Ce obtained. The values of C., and a 


are tabulated in table 2 for each I. The column titled 'closest to 
normal’ in table 2 shows the transformation that give values of ve and 


“a 


E closest to meeting the criteria set by equations 7 and 9. This shows 
that several transformations can be applied to the data but it would 
appear that the square-root is the most applicable. The square-root gives 


the value of Se and C, closest to normal whenever the square-root and another 


k 


“Aw 


transformation lie inside the limits. On a few occasions the values of C. 


and Cc. lie just outside the limits of equations 7 and 9. The limits 


seem to be between those defined by equations 11 and 12 which are not 


all that different for the normal condition. 
Op ro hs 0.3 (11) 
(y) -1.0<% < 1.0 (12) 


the values of a and Ch all fall within these new limits for the same 


transforming agent except for the month of July where CL. iststili toe 


high. Table 3 lists the month and the suitable transforming agent. 
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TABLE 3 


DETERMINED TRANSFORMING AGENTS FOR PRECIPITATION 


Month Transforming Agent 
JAN 0.5 
FEB 0.5 
MAR | 0.4 
APR 0.4 
MAY 0.4 
JUNE 0.4 
JULY O20 5 
AUG 0.5 
SEPT 0.5 
OCT 0.5 
NOV OAS) 
DEC 0.5 


YEAR 0.25 


aD 


Inflow Parameter 

Inflow to a lake sub-system can take one, or a combination, of 
the following basic forms; surface runoff, groundwater (which can appear 
above or below the lake surface) and streamflow, primarily from one 
lake to another. Surface runoff and groundwater components of the 
inflow parameter will be dealt with in later sections as they warrant 
separate treatment. 

The schematic representation of the lake system shown in Figure 5 
illustrates the path water would take if each lake were full and dis- 
charging its surplus water to the next lake in the series. The stage 
at which the lakes begin to discharge is controlled by the level of the 
overflow channel. This level can only be estimated roughly, at this 
moment, from topographical maps of each lake. These estimates are 
given in table 4. The water balance model should have a routine incor- 
porated into it to estimate, for a given time period, the amount of water 
a lake will contribute to the water balance of another lake. Losses 
will naturally occur in the channel connecting the lakes, particularly 
after a dry, no flow period. These losses can be attributed to the 
moisture deficit that occurs in any natural channel if it is left dry, 
to evaporation from the water surface and transpiration of phreatophytes 
and other plants growing in the channel. The practical estimation of 
these losses is impossible due to the number of affecting variables such 
as channel dimensions, materials composing the channel sides and bed, 
time of year, frequency and magnitude of discharges through the channel 
and the type and state of plant growth. Moreover, the magnitude of 


_these intermittent losses is negligible when compared to the inaccuracies 
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of measuring the other water balance components,in particular surface 
runoff. 

By considering the event that a lake has an excess to discharge 
it is justifiable to assume that the magnitude of water transferred 
from one lake to another is equal to the excess of the contributing 


lake calculated for that time period. 


TABLE 4 


MAXIMUM AND POSSIBLE MINIMUM LAKE LEVELS ABOVE ORDINANCE DATUM. 


(OBTAINED FROM 1:25000 TOPOGRAPHICAL MAPS). 


LAKE MAXIMUM LEVEL MINIMUM LEVEL 

Pls M ue hs M 
Miquelon 2515 766 2483 756 
Oliver 2505 763 2485 igo} 
Joseph 761.5 (poe) 
Ministik 2500 761.4 2480 ee, 
Cooking 2423 738 2400 foes 
Hastings 2414 735 Zoos 730 

Groundwater 


The flow of groundwater into and over the subsurface boundaries 
of the Cooking Lake Moraine in governed mainly by the topography of the 
area. To describe in detail the numerous groundwater flow patterns that 
could exist in such a complex situation as is to be found in the Moraine 
would require an extensive well monitoring program and is beyond the 


scope of this study. The Moraine is particularly complex due to the 
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diversity in surface and sub-surface materials and to the complicating 
aspects of permanent lakes on any groundwater system. The Moraine 

is a local topographical high and as such can be generalised to the 
typical prairie profile (Meyboom),1962, 1966), but it must be said 

that the local groundwater inflow and outflow of individual lakes is at 
the present time an unknown factor in the water balance of each lake. 

The prairie profile has been offered as a model of the typical 
prairie groundwater flow condition to which all observable groundwater 
phenomena can be related. The basic properties of the model are shown 
in Figure 7. Geologically the profile ideally consists of two layers, 
the uppermost being the least permeable, with a steady-state flow of 
water to a discharge area. The prairie profile can often be considered 
as a model for small scale systems, such as the typical knob and adjacent 
kettle, common to rolling prairie topographyg and as a model for wider, 
regional scale systems. The applicability of the profile has been sub- 
stantiated by numerous borings for both scale systems. 

It is often quite feasible to delineate between recharge and dis- 
charge areas because each has particular and distinctive characteristics, 
the most prominent being the change in potential with increasing depth. 
An increase in potential with depth indicates a discharge area and a 
decrease in potential with depth a recharge area. Flowing wells are 
found in discharge areas in many instances. 

When groundwater is discharged from the groundwater system to a 
surface system it is usually lost from the entire water system by 
evapotranspiration. Considerable attention has been given to evapo- 
transpiration and to the many surface features that are associated with 
it. The estimation of evapotranspiration, in particular, will be dis- 


cussed in the next Chapter. The remainder of this section will concen- 
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local flow system 


intermediate flow system 


(a) 


zone of depression beneath 
phreatophyte fringes 


shallow leakage 


Legend: Direction of groundwater movement man ete 


Watertable Rereie Be 


Boundary between flow systems of different magnitude - --~-~ 


FIGURE 8: Diagrams of Flow Conditions Near Permanent Lakes 
(after Meyboom, 1966) 
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trate on the surface features associated with discharge and recharge 
areas. 

The spatial investigation of various vegetal types, within a 
given watershed, will often yield valuable clues to the physical 
character of groundwater flow and will directly give a good indication 
of soil type (hence for runoff characteristics). The usual  obser- 
vations associated with discharge and recharge areas include the 
following, 

(a) the occurrence or non-occurrence of willow rings and the 

chemical character of water bodies centered within them 
when they do occur. 

(b) The distribution of halophytic communities with respect 

to the occurrence of local and regional flow systems. 

(c) The location of ponds and bogs with respect to groundwater 

flow. 
The occurrence of willow rings centered around higher areas, that is, 
recharge areas, indicates local discharge points within the recharge 
area, [Meyboom, 1966]. Each discharge point is characterised by a 
willow ring. Willows have a very low tolerance for saline soils, 
indicating their affinity for waters which have not travelled far in the 
groundwater system. The typical basin will have numerous willow rings 
in its' higher areas, wheras the lower regions are void of rings. The 
occurrence of willow rings in an indication of a local flow system 
superimposed on the regional groundwater flow system. [Meyboom, 1966]. 

Whenever there is a nett upward movement of mineralised ground- 
water, extensive saline soil areas occur. In many cases the transitional 


changes of groundwater and vegetation within a basin are from relatively 
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pure water and willow rings in the recharge areas to saline water and 
halophytic plants in the discharge areas. When the local flow system 
is not replenished by saline water formations, the groundwater is 
relatively fresh and saline soils fail to develop. However if there 
is still a net surplus of water the result will be the occupation of 
the area by fresh water phreatophytes. Several of the above features 
can be found in the Cooking Lake Moraine. 

In areas of hummocky terrain (see description of the surficial 
geology, Chapter 1) ephemeral water bodies have been found to function 
as recharge points during spring and early summer, and as discharge 
points during summer and autumn. Permanent lakes however have been 
shown to be areas of permanent groundwater discharge. [Meyboom, 1966, 
1967]. The four typical flow conditions near permanent lakes are shown 
in Figure 8 a-d. These four conditions are believed to be similar to 
the groundwater flow conditions found around the lakes in Cooking Lake 
Moraine. 

The first situation, figure 8a, is a spring condition where dis- 
charge from local and the intermediate flow system moves towards the 
lakes. Figure 8b is an approximate summer condition of seepage towards 
the phreatophyte fringe surrounding the lake. Insufficient recharge 
causes a deterioration of local flow and will result in shallow move- 
ment of groundwater from Lake A to Lake B (Figure 8c). Figure 8d 
shows the final deteriorated system for the autumn and winter condition. 

The numerous studies made of the groundwater flow condition near 
lakes demonstrates the dynamic nature of lake bodies and the groundwater 
flow condition. This implies the virtual impossibility of describing 


the nature of groundwater movement, in the vicinity of lakes, in terms 


Bh fii 
ee 


aoe ag Le 
axsids aaa caphikdne ibid ' awh Ye 
q we S Ris oe ah 
20 @obiaquasy oft wd LEEW: ditto oift sre! 4 L i 
sewidnei svods ofa to faievse |) me wo 
caibwrelt sdad’ i . 


OY, a i at 


; ir 
va 7 a, : 
Pritam ft 


tatoitiue sf? Yo wotsgkyoesb san) nieve go niee 
coljomt of bavel need sve ceaiatal “aw sia ea 
 gaukoq opradocs 8 B 
pe rommitye github sad 
AICI jmsctys] .aguadaekh weyevhaoggy: samatiiiog 20 ass ad of inode 
ewok ats adel snomemreq 289m ehokstbaod wolt \hdekar> 1902 ‘eit .[teet — 
oy wHilmiz sd-o9 Siam fed) ete enoktibnds spot seed? .b-s 8 ought ah 
sted gottoud af avdad an2 banons brat anotikbaes wol? wosubewong, Pre) oan 


conten 4) i 7 r 


-ath stedw oabethnos gnithe o at eO amight ,nobyautte sedtY air” v | 7 


Ssutatoeth es bas ees vites brs gobrqe 4 


aeod avad tevewol eodsl tasanoset =. nen 


ait ebtewot asvem agseye wold sancbonrise ma ef2 brie Iscof mort sgt 
abrawod sgsqgse lo poldtboos teuieve saadikorqge he af d8 emwglt waite” 
sysuiioos anstotiuent sale ect gabbauctive, syaket saydqoseo te edt 
-svem woitede at tiyesy iirw boas yor? Lapel to boksarateenab x ee 
bf sult C9 sxngi) & stad os A ldiiad woe) Secrmwbimmayg <4 a 
nokikbino wonky bos iMimtus sit eT mesage bedwrotra39b Ennk? od wh 5 3 
quad mobshhnes woll tesswbnuony ods te stem asthuss eyorsaun — 
is jevbrmeng sda) bas) sekbod sie to mane senso aes ven nc 
Saad in wat itdtewoqal feaovie sty esktgatt: enue’ eno 3b 


dial le nN 


jar 


| eal 
te ie 
ke - 


7 al: ' Lo ‘1 hi, 


ry. i oy mi . 4 : aan 


32 


of static analysis. In many cases it is safe to assume that there are 
sub-surface flow into the lake from the watershed and sub-surface flows 
between the lakes from the watershed. The quantities of water involved 
in each sub-surface flow system is believed to be negligible in compari- 
son with the surface runoff and evaporation components of the water 
balance of the area. A low quantity of groundwater flow between the 
lakes can be seen to be feasible if the geological cross-sections 
defined in Figure 2 and shown in Figure 3 are inspected. In Figure 3 

it can be seen that clay (an aquiclude) has extensively covered the 

area and overlies another low permeability material, shale. 

Springs, above and below the surface of Miquelon Lake have been 
reported [Nyland, 1970] and Duck's Unlimited mentions the presence of 
springs at Ministik Lake in the 1930's. It has been hypothesised that 
springs must be an important feature in the hydrology of Hastings 
Lake if this hydrology is to adequately explained. [Kerekes, 1965]. 
All of the above observations point towards the complexity of the 
situation to be found in the Moraine. Therefore to model groundwater 
flows properly the geologyof the area must be investigated in greater 
detail than is possible here. The groundwater condition is not a 
static situation in the case of the Cooking Lake Moraine and so a 
"lump' type system of modelling will not accurately reflect the stoch- 
astic variations within the groundwater system. Consequently the 


groundwater contributions to the water balance must be ignored. 


Surface Runoff 
In very general terms the water yield of a basin (which includes 


surface runoff and groundwater flow) is simply the difference between 
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the precipitation over the total watershed (land surface and open 
water surfaces) and the losses due to evaporation and transpiration. 
The combined loss due to evaporation from open water surfaces and 
the ground and transpiration by vegetation is an important component 
of the water balance equation which is difficult to quantify. 

Many techniques have been developed to estimate evapotranspiration. 
Some techniques rely on elaborate theories, such as the aerodynamic 
approach or eddy-correlation theory, others are based on these theories 
but have been empirically derived. Because of the importance of 
evaporation and transpiration in any water balance model a brief over- 


view of the main estimation techniques is given in Chapter 3. 


Summary of Chapter 


If the losses due to evaporation and evapotranspiration are 
known and the amount of precipitation that has fallen over the water- 
shed is measured, for the same time period, it is possible to give a value 
for the expected water yield. The spatial variation of precipitation, 
in particular, over any area is one of the most variable parameters 
in the water balance of that area. As the spatial measurement of 
precipitation over any appreciable area is inpracticable it is necessary 
to adopt point measurements as a representation of spatial precipitation 
amounts. The problems encountered and the errors associated with adopt- 
ing point measurements as a model of regional precipitation has been 
the topic of many papers and discussions. [Ferguson and St rr in 1967]. 
Groundwater contributions to the water balance of the Cooking 
Lake Moraine has been shown to be a dynamic problem with a considerable 


variation between the seasons. These contributions however, are not 
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as important as the surface runoff or precipitation components of the 
water balance so the assumption of negligible groundwater input or 
output from the water system is acceptable. 

Like precipitation, evapotranspiration estimation is limited by 
the extend and practicability of the variables measured, the type 
of variable measured and the reliability of the measuring techniques. 
In the long run the estimation technique adopted to give a measure of 
surface runoff is determined by the availability of information. 
Information that is available consists of records of precipitation 
(total daily at Edmonton), mean daily temperature (maximum, and minimum) 
for a relatively long period, that is, greater than forty years. 
Information such as maximum and minimum lake levels are obviously 


limited in accuracy because of the means by which they were obtained. 
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CHAPTER 3 


ESTIMATION OF EVAPORATION AND EVAPOTRANSPTRATION 


This Chapter deals with the estimation of open water evaporation 
and evapotranspiration from land surfaces. It is necessary to deal 
with this topic in the detail given here because of its" relative 
importance in the hydrology of any area and it is especially important 
when investigating a lake district such as the Cooking Lake Moraine. 
Some of the methods outlined are only suitable for the more elaborate 
or extensive water balance studies. Consequently these methods are 
only of academic interest and have been included in this dissertation 
only for completeness. 

The selection of a particular method or empirical formula is 
based primarily on the degree to which the area under study has been 
monitored. Whether or not the various independent variables of each 
method have been measured and are available will determine the selection 
of the technique incorporated into the model to measure the evaporation 
or evapotranspiration component. 

There are several methods of estimating evaporation and trans- 
piration, the accuracy of each method varying considerably from one 
location to the next and with respect to eachother. The more commonly 
used methods include [Munn, 1961]. 

(a) The Water budget (or storage equation) 

(b) The Energy budget 

(c) The Eddy-Correlation method 

(d) The Aerodynamic-profile method (humidity and wind 


velocity gradients). 
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Combinations of the above four methods form the basis of many 
of the empirical formulae which have been derived and are in use today. 
The more popular of these formula are those which have been derived 
by the following experts and the formula's usually carry their name 
(1) Penman; (2)... Cherntnwaite and (3) Ture 
Many other formula have been derived by combining properties of the 
above estimation methods. The more recent formula's tend to improve 
the results when compared with the results of classical formulas but 
are usually only applicable to the conditions found at one locality. 
This type of local formula include those of, Blaney and Criddle, McIlroy, 


Rohwer and Meyer and Tichomirov. 


Methods Used to Estimate Evaporation 

There are several formal approaches used in the estimation of 
evaporation and transpiration, two of which are equally applicable to 
evaporation from a free water surface and evapotranspiration from 
vegetation. These two are the Water-budget method and the Eddy-correlation 


method. 


Water-Budget Method 


This method generally involves the use of a storage equation of 


the form given by equation (14) 


P+IT+u=E+0O+A8S (14) 


where 


P = total precipitation falling on the area 
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I = volume of inflow 
u = volume of underground inflow or outflow 
E = evaporation, or evapotranspiration from the surface 


of the area 
0 = outflow volume 
AS = change in storage (in reservoir or groundwater) 

The water-budget method can be applied both to evaporation and 
evapotranspiration but the estimation of the latter is more difficult 
than that of the former. The main reason for this is the difficulty 
in estimating the change in storage, particularly groundwater storage 
changes, and losses to deep percolation. Changes in water-table level 
can be easily monitored but translation of such fluctuations into an 
actual change in storage is difficult. 

The units of equation (14) are usually inches or millimetres over 
a given area for some convenient time interval. One main disadvantage 
of this method is that each variable can only be measured to a particular 
degree of accuracy and accumulation of these errors in equation (14) 
would result in poor estimates of evaporation. To illustrate this point 


considertthe F.D. Roosevelt Lake, formed by the Grand Coulee Dam, on 


the Columbia River. The annual evaporation of this lake has been estimated 


at forty-two inches, or equivalent to 290,000 acre-feet. If the reser- 
voir outflow (estimated at about 70 million acre-feet) was measured 
EO an, eccuracy GL 57, Lue eCrror would be apo + 3.5 million “acre-teer 


or more than ten times the total evaporation. It can be seen that for 


certain large catchments or large lakes, an accurate water-budget approach 


to estimating evaporation or evapotranspiration is not practical. Many 


safeguards are necessary to ensure that all variables are measured to a 
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realistic degree of accuracy. If this accuracy could be achieved the 
water-budget method would be ideally suitable for evaporation estimation 


particularly for the smaller watersheds. 


Energy-Budget Method: 


The energy-budget or heat-budget was first used by Schmidt in 
1915 to compute evaporation from oceans. Many other investigators have 
attempted to use the technique to compute evaporation from bodies of 
water of all sizes and to estimate evapotranspiration. Only in the 
last decade has the necessary instrumentation become available to 
permit the measurment of certain variables in the technique. 

The energy-budget method is based on the principle of the conser- 
vation of heat-energy within a body of water. This conservation principle 
states that a balance must exist between; (a) insolation, (b) heat 
transferred from the water surface by radiation, convection and conduction, 
(c) heat energy acquired or lost in raising or lowering the temperature 
of the water, and (d) heat dissipated or acquired by evaporation or 
condensation. If items (a), (b) and (c) can be measured, evaporation or 
condensation may be estimated. 

The water-budget was used by Anderson in 1950-1951 as a control to 
test the energy-budget method for estimating evaporation from Lake 
Hefner, Oklahoma (Anderson, E.R., 1954). He concluded that the method 
gave satisfactory results for periods of ten days or more. Unacceptable 
results were obtained for periods less than 10 days for Lake Hefner. 
Generally periods of one or two weeks have been used for most energy- 


budget estimates. 
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The use of the energy-budget equation to estimate evapotranspiration 
is not applicable because the evaporating or transpiring surface is 
three dimensional. For evapotranspiration measurements using an energy- 
budget appraoch a more complex energy relationship is required. Assump- 
tions to the extent that energy losses due to precipitation, Snowmelt 
and runoff are negligible and that energy storage is minimal, are 
required if a workable equation is to be found. The type of equation 


derived is similar to equation 15. 


RY aS 
outa of, 
where: 
E = evaporation 
RO = net radiation flux 
S = soil heat flux 
8 = Bowen's ratio, which can be defined as the ratio of 


the loss of upward energy flux as sensible heat to 
the energy flux used in evaporation. 
To reach reliable conclusions from energy-budget studies, detailed 
measurements under different weather conditions are required, preferably 
continuous records of the various variables in the energy balance 


relationship. 


Eddy-Correlation Method: 
This method, also known as the eddy-flux or eddy-transfer method, 


is equally applicable to evaporation from a free water surface or to 
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evapotranspiration.' It employs the measurement, at some arbitary 
point, of vertical turbulent fluxes in the atmosphere of both velocity 
and water vapour. The basic principles of eddy-correlation were first 
recognised in 1951 (Swinbank, 1951). 

From a theoretical point of view, this method is quite straight- 
forward. However, the measurement of simultaneous fluctuations of 
the deviations from the mean of the vertical wind velocity and specific 
humidity, requires an instrument sensitive enough to measure these 
rapid changes and,sturdy and stable enough for continuous use. At 
the present time no instrument has been developed that is capable of 
measuring only the eddy-flux. Instruments do exist however that do 
measure total-flux and eddy-flux combined with the flux due to mean 


motion. 


Aerodynamic-Profile Approach: 

The two premises below are the base of this method [Sutton, 1953] 

(a) If a moisture gradient exists in the air, water vapour will move 
towards points of lower moisture content. 

(b) The rate of movement of the water vapour is accentuated by the 
intensity of turbulence in the air. 

The profile technique concerns itself with the turbulent transfer 
of water-vapour between two levels in the air at a small distance above 
the evaporating surface. It is usually assummed that: 

(1) Diffusion in the forward direction is negligible in 
comparison with forward transportation by the wind. 


(2) Steady-state conditions exist, and 
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(3) There is.an infinite line source of water, that is, the 
problem is two-dimensinvnal. 

An equation can be derived by analogy with molecular and other well 
known transfer processes. Accepting assumptions (1), (2) and (3) are 
valid and that there are spherical eddies with a logarithmic wind 
profile above the evaporating surface, an equation of the form given 
by equation 16 can be derived. This is in effect the classical aero- 
dynamic-profile equation derived by Thornthwaiteand Holtzman. (Thornth- 
waite and Holtzman, 1939) 

12 
e(g, - 8,)(u, - u,) 


[log (Z,/Z,)] 


where 
E = evaporation in cm/s 
e = density of air in neem: 
85 and 8, = specific humidities at heights Zo and 2, respectively 
u, and u, = wind speeds (cm/s) at heights Z5 and zy respectively. 
and K = Karman's constant = 0.41. 


It should be again noted that this equation is strickly valid for 
neutral (steady-state) conditions. Under other conditions, such as 
periods of intense evaporation, use of the equation will give erroneous 
values due to the breaking of the logarithmic wind profile assumption. 

Special equipment is required to make the necessary measurements of 
the variables in equation (16). Very accurate sensors have to be used 
to measure the small gradients of wind speed and vapour pressure over 


the short height interval. The temperature sensors must be able to read 
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correct to 0.05°C and humidity sensors to 0.01 mb. The anemometer 
must have minimum internal friction and a minimum time lag in response 
to wind fluctuations, they should (ideally) have a starting speed of 
less than 10 cm/s. 

The Thornthwaite-Holtzman relationship is applicable to aerodynami- 
cally smooth surfaces such as lakes and reservoirs where the level of 
effective evaporation is taken as the water surface itself. On land, 
if the vegetation height is much different from the heights above the 
soil surface at which vapour pressure and wind speed measurements are 
made, the wind profile requirement can only be satisfied by applying a 
"Zero displacement correction". This is due to the boundary from which 
turbulent transport is effective and is at some distance above the soil 
surface. An adoption (Rider, 1954) has been made to Thornthwaite's- 
Holtzman's relation to include a term for zero displacement. (Rider, 


1954). 


Empirical Formula: 

Many attempts have been made to produce satisfactory formula for 
the estimation of evaporation. Usually they try to estimate evaporation 
from open water-surfaces but some others estimated potential and actual 
evaporation and transpiration from vegetated surfaces. Most of these 
empirical formulae assume a constant supply of water hence estimates 
potential evapotranspiration. 

The formula developed usually use data from meteorological measure- 
ments to compute evaporation or transpiration and are based on a 
combination of two or more of the estimation methods described above. 


Perhaps the most widely used approach for estimating evaporation from 
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meteorological factors is that based on a combination of the aero- 
dynamic profile and the energy-budget methods. Of this type the 


Penman equation is the best known. 


Penman Equation: 

(Penman, 1948, 1952, 1956). Penman originally published his 
formula in 1948. It was used to estimate evaporation from free water 
surfaces from meteorlogical data. The formula requires the assumption 
of unlimited availability of water for evapotranspiration. Review of 
the technique in 1952 and 1956 resulted in lower and more accurate 
estimates for evaporation. Penman's approach is based on the measure- 
ment of the following four meteorological variables. 

(1) Duration of bright sunshine 

(2) Air temperature 

(3) VAL humidity 
and 

(4) Wind speed. 
These four factors reflect most of the terms in the energy relationship 
(equation 15) but they do not take into account the energy reaching 
the body from the ground, rivers, precipitation and snowmelt for example 
nor do they show the effect of energy stored in the waterbody. 


The simplified version of Penman's equation is given by (17) 


cS H + E 
E = oo (17) 
miipor 
i ) 
where 
E = evaporation from surface. (E) for open water and Ee 


for evapotranspiration). 
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A = slope of standard vapour and pressure curve at mean 


air temperature T 


me 
yY = the constant of the wet and dry psychrometer equation 
H = net radiation reaching the hody 

ES = an expression for the drying power of the air 

x = an expression indicating surface type. For open water, 


x = 1. For vegetated surface x > l. 


For open water, therefore, equation (17) becomes 


(H+E) 
ae Ne (18) 
(~-+1 
. ) 
Penman introduced a function which would vary with time of year and 
would enable direct estimation of evapotranspiration from estimated 
open water evaporation. This is given by equation (19) and f varies 
irom 0.6 to. 0.8 
a =fx ES (19) 


It is not always possible to measure terms such as net radiation 
so Penman and others have set out empirical expressions for these cases. 

E. in equation (18) is an expression relating the wind speed at a 
height of two metres above the evaporating surface and the saturation 
deficit (e, - e, of the air at this height. Penman's formula for 


potential evapotranspiration is given by equation (20) 
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(— - H+ E,,_) 
AT 
a oe 
—+ 1] 
CS ) 
where 
S (igs we, 
EAT = 0e3o° CLF 100? 6&4 Qe 


Thornthwaite Equation; ( Thornthwaite, 1939, 1948, 1954) 

The equation was developed in the United States using a statistical 
study of available observations. The dominant parameters in the appraoch 
are temperature and length of day. Thornthwaite intended that temperature 
would be a reasonable parameter to integrate the balance of radiation 
used for heating the soil and air and radiation used for evaporation. 

The Thornthwaite model is based on calculations of potential evaporation 
or transpiration, then, on the basis of a series of assumptions, the 
estimates of monthly runoff are made using established empirical rules. 

Equation 21 lists the general Thornthwaite Equation for estimating 


potential evapotranspiration. 


E = C, Ti C21) 
where 
E = evaporation or potential evapotranspiration (water un- 
limited), 
C = coefficient 
T = monthly mean temperature ‘Gur and 
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The constants a and c both depend on the location of the area. 


The exponent 'a' can be evaluated in terms of the annual heat index, 


ir, as 
Dp =o -6 _2 
aa=- 67.5 x 10 Li Ge LG Lie OL9T + 05492 (22) 
in which 
iZ TF d ips Wi 
I =} eae (23) 


m=1 


Thomthwaite simplified this relationship for the case of 12 hours 
of sunshine per day for a 30 day months. Equation (21) reduces to, 


i i 9 ed ite 


E= 1.62 [ ; 2h) (24) 


For all other cases, tables are available which relate a correction 
factor for daylight hours to latitude. This correction factor is the 
constant 'c' in equation (21). Although the Thornthwaite equation is 


somewhat cumbersome to manually operate it does reduce to a simple nomo- 


graph [Gray, 1970] and lends itself easily for computer application. 


Blaney and Criddle Method: [Blaney and Criddle, 1950] 

This is similar to the Thornthwaite method in that it assumes that the 
heat budget is shared in fixed proportion between heating the air and 
evaporati'n. The formula was developed for the arid western United 
States and estimates comsumption use based on measurements of temperature 
and hours of sunshine. The monthly consumptive use, cu, is found by 
multiplying the mean monthly temperature, TT. with the monthly percent of 


annual daytime hours, p, and a monthly crop coefficient, k. 
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(Equation (25)) 


k tT. See 
a 8 00 eo (25) 
tee 
where anes monthly consumptive use factor 


The seasonal consumptive use is obtained by summation of the 
relevant monthly values. The coefficient, k, is an average seasonal 
consumptive use coefficient selected for different crops. 

The last three formula are based on the assumption that the rate of 
water use is unlimited, a condition seldom realised in nature. Several 
drawbacks exist for each type of method and these should be noted 
(Gray, 1970). 

(1) The Thornthwaite and the Blaney and Criddle methods use 
the same data to arrive at an estimate of annual or 
seasonal consumptive use. Generally they both tend to 
give an overestimation of water use in the early growing 
season and an underestimation in mid-season, unless an 
appropriate crop factor is applied. 

(2) The Blaney-Criddle method gives more reliable estimates 
of seasonal use than the Thornthwaite method for arid 
regions. They both are not as realiable as Penman's 
method for humid areas because of the absence of a humidity 
term. 

(3) The Thomthwaite and Blaney-Criddle methods should not be 


used for estimations of water use for short periods of 
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time as no allowance is made in the methods for 
variations in wind and relative humidity. 

(4) The Penman equation is hindered in its' application due 
to the number of meteorological observations required, 
however, it does provide the best means of estimating 
potential evapotranspiration. Some care should be 
used when applying to arid regions as advected energy may 


by significant. 


Turc Equation: [Turc, 1954, 1955] 

This formula was developed by assuming that water supply is not 
unlimited, but it is still basically a temperature formula. Turc 
carried out a comprehensive survey of general experience and found that 
usually evaporation is greater in wetter years than in drier. The 
first formula was designed to relate annual evaporation from catchment 
areas to precipitation and air temperature, the latter giving an indication 


of evaporation opportunity. The first Turc formula is given by equation 


(26) 
ee +52 (26) 
[0.9 + G) ] 
where 
E = evaporation (mm/annum) 
P = precipitation (mm) 
L = 300 + 25T + 0.05T° 
T = mean air temperature (°C) 


A more complex formula was latter developed for shorter time 
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periods were varying levels in soil moisture and crop factors are more 


pronounced. (Equation (27)) 


Peta + V 
E = i Pee eee (27) 
+ ——— — 
where 
E = evaporation in a 10-day period (mm) 
P = precipitation in the same period (mm) 
a = estimated evaporation from bare soil during the same 
period assumming zero precipitation (a < 10) 
V = a crop factor 
2 = evaporation capacity of the air 
Dore e % 
x ( ) si 
16 
eae : ; 2 
R., = incoming radiant energy (cal/cm ) 
and al = mean air temperature over 10-day period (oe) 


The Cooking Lake Moraine has limited data on hydrological and 
meteorological variables. The results of this fact is the necessity 
of adopting procedures and techniques which are not the most accurate. 
Only two meteorological variables are available for any considerable 
length of time at a station close enough to be relevant to this study. 
These are the ninety-one years of record of maximum and minimum daily 
temperature and total daily precipitation recorded at Edmonton. This 
paucity of suitable records of meteorological variables limits the 
evaporation estimation technique to one which needs only temperature 
and precipitation. There are several such methods available but that 


devised by Thornthwaitehas been selected. 
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CHAPTER 4 


MONTE CARLO METHODS 


Introduction 

The Monte Carlo method can be defined as a technique which consists 
of generating a long series based on the characteristics of an available 
Sample series. The method is also known as the experimental statistical 
technique. The definite process, in the Monte Carlo Method, is replaced 
by a random or stochastic process which will give the same result. It 
is worthy to observe that although a longer period of generated values 
are available, the information contained within the series is no more 
or no less than that of the sample series. Even with this limitation 
the technique has found various degrees of success in solving the numerous 
hydrologic problems [Sudler, 1927. Swanidze, 1964. Yevjevich, 1965]. 

It has of course had a longer and wider use in solving other probability 
problems such as, nuclear disintegration, cosmic clouds, telephone or 
traffic studies and in making random choices in some decision making 
process, say in playing a game. It is in fact probabiy correct to say 
that the experimental statistical method is as old as probability theory 
itself. 

The first recorded example of the Monte Carlo process was probably 
that of Count Buffon's needle experiment of 1773. In this experiment 
Buffon observed that if a needle of length less than or equal to unity 
were tossed at random onto a horizontal surface ruled with equally 
spaced lines (at unit spacing) then the probability of the needle cross- 
ing the line is equal to 2L/17, where L = spacing. Buffon therefore 


reasoned that he could determine the value of 1 experimentally by making 
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repeated trials. Another classic example of the Monte Carlo process, due 
this time to Roger Pinkham, is the experimental estimation of e. This 
experiment is based on the observation that if 2k numbers Xs» Pee"). Zk, 
are drawn in sequence from a random, equilikely source, then the pro- 
bability that they are all in ascending order, that is, Xs Xo» Xqe see 


x is equal to 1/(2k)! If the probability of a failure, of an ascending 


2k 


sequence, on the odd number 2k + 1 is considered, that is, the difference 
id OAS) = 1/(2k+1)! 


it is possible to obtain the total probability that a sequency of draw- 
ings of random numbers from an equilikely source will produce a rising 
sequence that ends with an even number of numbers. This probability is, 


co 


cn iy! ~ Girt) ae Cai! ~ Gey] * ai 
From an experiment consisting of 252 runs a value of l/e of 0.381 
was obtained, this represents an error of 3.5% in the estimation of e. 
Student in 1908 developed the idea of random sampling to estimate 
given distribution functions. Decks of cards and other sampling devices 
have been used to generate non-historic river flow patterns (Sudler, 
1927) to achieve a mass diagram analysis to develope probability distri- 
butions of reservoir capacity. This was probably the first serious 


application of the Monte Carlo process to solve a hydrological problem. 
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Sources of Random Numbers 

The first source of random numbers that likely springs to mind 
is some natural, physical phenomenon such as, cosmic rays, nuclear 
disintegration or perhaps sampling at random an alternating frequency. 
This technique of using a physical source has found many applications 
in past studies but has generally been superseded by computer methods 
which use some form of recurrence relation. A recurrence relation 
is where each successive number in a series is formed from the pre- 
ceding number or in some relationships from several of the preceding 
numbers. Usually some form of a algorithm of arithmetic operations 
is applied. The most important advantages of the recurrent technique 
over the physical sampling are the repeatability of the experiment 
and the non-reliance on any physical properties particularly the 
stability. 

The numbers produced from a:scource number via a recurrence relation 
are predictable. In the true mathematical sense of random, the numbers 
Simulated by the recurrence relation are not pure random numbers. To 
delineate then between the impossible, pure random numbers and the 
predictable random number, the term psuedo-random is used. The numbers 
therefore generated by a recurrence relationship are psuedo-random num- 
bers. In practice it is of no consequence if the numbers are predictable, 
what is of importance is whetheror not the use of the numbers is 
random. 

Von Neuman has developed an algorithm method to generate psuedo- 
random numbers based on selecting the middle digits of products (Von 


Neuman, 1951). This particular generator is an example of a general 
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class of psuedo-random number generators which are most often used 
to-day and are commonly called multiplicative congruence generators. 


These are defined by equation (29) 
ny k 
Xx = x o (MOD 2°) (29) 


or as "X41 is congruent to x modulo okn which means the difference 


between X and X is divisible by 2k. 
n+ n 


1 


Selection of p and Starting Value 


First, it is obvious that if the starting value (X,) is even, that 
is, divisible by 2, all the products will have at least one zero at the 
end. This effective wastes one bit of the machine capacity. Consequently 
Xx. should be odd. If p is even,zeros will accumulate on the right hand 
Side until 


A ee k 
xy fey X, = 0 tod (30) 


which results in all numbers past this point being zero. Therefore o 
is also an odd number. 

Odd numbers can be written in one of the following four ways: 
(i) GOte e330 (2) “8t + 13... 03)" bte- aewende) (St. =—3.. for some 
interger t.The value suitable for the integer t has been widely inves- 
tigated (Hamming 1971, 1973) but would not be productive to reproduce 
the theory behind selection of t here. 

All of the multiplicative congruent generators produce a series of 


uniformly distributed numbers. As this distribution is not one of the 
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more useful hydrological distribution the problem arises of how to 
obtain another distribution from a sample uniformly distributed. In 
principle the following device works. If the wanted cumulative dis- 
tribution F(Y) of the distribution f(y) is equated to the cumulative 
distribution X of the flat or uniform distribution, solution of the 
equation will give numbers with the required distribution function. 
This is shown by the following equations. 
m8 


W 
ie dx =X al Ue ee ee) (31) 
(e) Oo 


Applying the inverse operator ria gives 


Fo F(y)} = y = Fo (x) (32) 


As an example suppose the exponential distribution is wanted, that is, 


2 ? ethen since, 


f(y) =e” 
4 iy) ay 
F(y) = f e dy sick Ay ¢ 
Oo 
As, F(y) = X by definition 
ey =] - xX 
or y = -ln (1 - x) 


By using the sequence of psuedo-random numbers generated by a suitable 
generator, Xs the series ys = In (x, ) is obtained, with Ys having an 
exponential distribution. By inverting the cumulative function the 


normal distribution can be obtained except that the inverse function has 
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to be approximated, usually by auxiliary tables. An alternative method, 
used in this study, is to call on, the central, limit theorem. In this 
method several psuedo-random numbers are added from a uniform (or any 
other) distribution. Since each number supplied from the uniform distri- 
bution is independent of all other numbers generated and the variance 


of the series is, 


Ap 
at =f «- =) dx 
fe) 
(x - 1/2) 
3 O 
Z 
oy a2 (33) 


Therefore it can seem that the sum of twelve distributions will have 
a variance of one. Thus the common rule is to add twelve numbers from 
the generator and subtract six from the sum to get a normal distribution 
with mean of zero and a variance of unity. 

A subroutine called grand is available in the University computer 
system, specifically the IBM Scientific Subroutine Package. The 
program has been writeen in assembly language but is essentually the same 
as subroutine RAND listed in appendix 2. This subroutine (grand) was 
used to generate a given length of precipitation and temperature data 
which follows a normal distribution with supplied mean and standard 
deviation. Grand computes twelve uniformly distributed numbers by the 
power residue method and with the central limit theorem obtains a psuedo- 
random number x conforming to the normal distribution. The random 


number x, will have the same mean and standard deviation as that supplied 
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to Grand. It was noted earlier that the starting value, Xo should be 
odd if machine space is not to be wasted. In fact any odd number 
between 1 and 2**(31-1) can be used for the starting number in the 
University of Alberta machine. Any number can therefore be applied 
and facilities do exist to enable the programmer to add this but sub- 
routine Grand supplies the starting value of 524287 if no other odd 
number is specified. This number was used for generating samples in 


this thesis. 
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CHAPTER 5 


LAKE LEVEL ANALYSIS 


As the title of this Chapter indicates it describes the analysis 
of lake level changes, obtained by the water-balance model, due to 
the variability of meteorological parameters in particular. The 
Chapter is made up of two parts. Part one deals with the proving of 
the model by estimating lake levels from the historical records of 
temperature and precipitation, and comparing these with lake levels 
measured. Part two extends the analysis by using data generated via 
the Monte Carlo Method (Chapter 4) and shows the possible large fluctua- 
tions that occur due to the natural random fluctuations in precipitation 


and temperature. 


Analysis Using Historical Data 


The method devised by Thornthwaite has been selected as the tech- 
nique for estimating evaporation and evapotranspiration. Laycock has 
shown the applicability of the Thornthwaite technique for numerous water- 
balance studies within the Canadian Prairies. (Laycock, 1964, 1967, 1971, 
1973). In his 1973 paper to the Symposium on the lakes of Western 
Canada, Laycock reviewed the basic water-balance characteristics of the 
Cooking Lake Moraine and used Thornthwaites formula to estimate surplus 
and deficit amounts for Edmonton data. 

Thornthwaite's equation is based on the mean monthly temperature. 
This gives an estimate of monthly potential evapotranspiration given an 


annual heat index composed of monthly heat indices. With known monthly 
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precipitation, storage capacity of the ground and antecedent storage 
condition estimates of monthly moisture surplus or deficit can be com- 
puted. Surplus conditions only occur whenever the difference between 
precipitation and evapotranspiration is positive. The surplus obtained 
can either runoff entirely or in part depending upon the antecedent 
storage condition. The runoff criteria can be defined by the following 


set of equations 


Sep Ste di Binaeks (34) 
S - SC S nee SO (35a) 
n n 
R = 
0 s, < 8¢ (35b) 
where 
Sn = storage condition of the soil horizons during month 


n, n=l, 12 (cms) 
P = total monthly precipitation for month n (cms) 
E = Thornthwaites estimate of monthly evapotranspiration (cms) 
SC = storage capacity of the ground (cms) 


R = runoff (cms) 


Le SF > SC there is runoff equal to the difference between the 
amount of water needed to fill the storage capacity and the amount of 
surplus. If runoff is greater than zero the storage condition of that 
month is set at SC. If Sa < SC there has been some depletion in the 


storage condition of the ground, that is, evapotranspiration has exceeded 
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precipitation. As each month's runoff or depletion is based partly 
on the previous month's storage condition, the Thornthwaite procedure 
gives an on-going method which models dry and wet periods. 

Master (appendix 2) is a computer program, written in Fortran, 
which produces estimates of monthly evapotranspiration given records of 
monthly mean temperature and total monthly precipitation and starting 
conditions. Various storage capacities can be assigned to the different 
geological, vegetal or topographical groups found in the Moraine. To 
accurately map the different storage conditions that exist throughout 
the Moraine would require a detailed investigation beyond the scope of this 
study. Consequently a typical or average value of four inches has been 
adopted to represent the unified storage capacity of the watershed. 
Laycock (Laycock, 1971, 1973) has shown that a storage capacity of 
four inches is representative of the majority of the watersheds in the 
Prairies and has applied this capacity to the Moraine area. A value 
of ten centimetres was adopted for this study to conform with S.I. 
units. A record of surplus and deficits for storage conditions (in steps 
of two centimetres) from 2 to 10 cms, and for 15 and 20 cms has been 
produced and are available if they should prove to be applicable later 
or when more input becomes available. 

The analysis of the fluctuations of lake levels using historical 
data as input was first attempted by assumming a constant lake area to 
contributing basin area ratio. The results of this analysis showed that 
it was not correct to assume a constant ratio. By an inspection of the 
area-elevation curves (appendix 1) it can be seen that the value of the 


ratio of areas is not constant and can vary considerably. Any sudden 
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increases in lake levels could therefore not be accurately modeled 
if the area ratio were set at a constant value. This ratio, D, is 
defined as the ratio of lake surface area to the effective drainage 
area. Ratio D was changed by unit increments from a low of zero 

(no runoff at all) to 5, believed to be the maximum value. Figure 9 
shows the lake levels of Cooking Lake for the six values of D and 

an arbitrary starting value of 100 m, also assummed to be the maximum 
lake level. 

Interpretation of Figure 9 shows the fluctuations that have 
occurred due to water deficiency and surplus over the period of record. 
This rapid rise in lake level, reported during 1900 and 1901, is graphi- 
cally shown. This sudden rise was due to the very heavy precipitation 
recorded during two consecutive years. The general "good condition" 
reported from the turn of the century to 1930 is shown by fairly stable 
lake levels. From 1930 to about 1970 a general decline in lake level 
can be seen for the higher values of D. It no runoff occurs, that is, D 
equals values of D zero, Cooking Lake would dry up completely in about 
forty years, assumming it starts from a full state Several times the 
lake reached a maximum level for values of D greater than four. From 
this analysis the author expects that the actual ratio between lake 
surface area and effective drainage area is around three to four. 

The adoption of a constant value for D is acceptable for a general 
appraisal of the lake level situation but it is not valid for the physical 
case. As the lake level increases and decreases over the months and years 
it is obvious to assume that the value of D will respectively increase or 
‘decrease. The next logical step is to try and model mathematically this 


change in D. 


Area elevation curves are given in Appendix 1 for the major lakes. 
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in the Cooking Lake Moraine area. In each case the area~elevation curve 
has been approximated by one or as two straight lines depending upon 
the individual case. The equation of a typical line is given by 


equation (36) 


SA =A.C +B (36) 
where 
SA = surface area of the lake in acres for a given surface 
elevation, C, which is in metres above datum. 
A and B- = constants defining the linear relationship between 


SA and C. 
The contributing area, CA, can be thoughtof as the effective drainage 

area and defined as that area of the total watershed which will contri- 

bute to runoff. This can easily be shown to be equal to the difference 
between the total drainage basin area(TDB) and the lake surface area (SA) 
multiplied by a factor (F) which depends upon the areal extent of depressions, 
marshes, phreatophytes etc, that is, those areas which do not contribute 

to runoff except possibly in very wet years. Mathematically the contri- 


buting area can be defined as, 


CA = (TDB - SA) x F (37) 


Table 5 lists the drainage basin areas and the linear relationship 


between surface water area and lake level 
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TABLE 5 
LINEAR RELATIONSHIPS BETWEEN LAKE SURFACE AREA (SA) 


AND LAKE SURFACE ELEVATION (C) 


TOTAL DRAINAGE 


LAKE BASIN AREA RELATIONSHIP RESTRAINT 
TBD (ACRES) 
Miquelon 14660. SA=459./760*C-348607 C<762.8 
#1 and #2 
SA=270.848*C-204490 C>762.8 
Oliver 9600 SA=361.13*C-272168 C</56.28 
SA=270.85*C-203895 C>756.1 
Joseph 8400 SA-317.30*C-238981 RELA 
Ministik 22360 SA=1128.53*C-850912 ALI*C 
Cooking 46080 SA=3830*C-2807267 C27 3562 
SA=1459*C-1064344 C>735.2 
Hastings 26240 SA=361.13*C-263440 Gs (6552 


SA=816.63*C-598305 C>735.2 
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elevation. The restraint listed in column 4 of Table 5 is the elevation 
at which the area-elevation curve changes slope thus requiring a 
different straight line definition. The lines were all fitted by eye. 
The determined contributing area will, for a given lake surface 
area, provide D times the surplus runoff estimated by Thornthwaite's 


procedure. That is, for R > 0 


ALL = D x (R) - PE (38a) 
Be. (TDB - SA) if 
re oye ammaene Jo Fx R- PE (38b) 
and 
ALL = -E for< Rk <.0 (39) 
where 


R =run off 


ALL = change in lake level when surplus runoff occurs 


PE estimate of potential evaporation 
and all other terms as defined earlier. 

It is worthly to note that assumed maximum and minimum lake levels 
have already been given in Table 4. For completeness, the relative 
areas of forests and total watershed area for each drainage basin is 
given in Table 6. 

For each major lake in the system a computer run (using MASTER) 
was made and levels for each month of the year estimated. The levels 
in October for each lake have been plotted with help from the Comput- 
ing Services subroutine, CGPL. The plot of estimated lake levels using 


historical data for each lake are similar in shape but the fluctua- 


tions can be larger or smaller depending upon the lake in question. 
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Actual records of lake levels are available for Cooking Lake and the 
plots of lake levels for this lake have been analysised. The plots 

of the other lakes are given in appendix 4, and for Cooking Lake in 
Figure 10. The other variable in the determination of lake level changes 
is the factor, F, which measures the extent of contributing area 
relative to total area. The factor takes into account the losses to 
runoff due to depressional storage, snow drifting and phreatophytic use. 
Again Master was used to estimate lake levels with F = 0.25 up to 

F = 1.0 incrementing in steps of 0.25. Figure 10 shows estimated 

lake levels obtained with F = 0.65. Figure 11 shows the October lake 
levels estimated by using the four values of factor 'F'. It can be 
seen that the recorded lake levels seem to fall between the estimated 
lake levels when F = 0.5 and F = 0.75, hence the selection of 0.65 

for 'F' when producing the levels shown in Figure 10. 

Records of actual lake levels from 1956 to the present time are 
available for Cooking Lake. As a value for the lake level in 1885 is 
unavailable, it is necessary to assume the level and to relate the esti- 
mate with the measurements for thet particular starting value. After 
several trials a starting lake level of 736.5 m above ordinance datum 
was found to give the closest approximation when comparisons were made. 

The fluctuations, illustrated by Figure 10 can occur sharply, 
particularly sudden rises in level. These will occur whenever an excep- 
tionally wet year or a series of above average rainy years are encountered. 
The estimated values of the Cooking Lake level follow the measured 
fluctuations quite reasonably, but do debiate on a few occasions consider- 
ably from the actual lake levels. More records and surveys would both 


improve the estimates and give a better test of the model. The devia~ 


a Ui etn sR ba: 


} a ' ' ith 
sai? has ates 5 gnidoos 


; ¥ : + * : z 

siolq sat bog teyinad wood sve 

—_ 7 : a U 

bt siot’ gatiieoD sat tie .+ xtbnseqee see: 

| ; X 

asgnado Level sist Yo notieatarsdab ons ne, “on 


cle ° 
2 


yw 
g4ue gataudieinoo to Wweiwxs 95 ers 3890 


' Ly 7 i ao 
ag@eol #79 Gatinoss ‘otnt eeuel yo4SB1 Set 
» i 


;o 2 f See TD wou 
ay ¢ jtw elsvel sos 
; E +e fi ; ; p 
‘i + ta 3 eS 


yFamlyes sb mwtad List i mesc A ‘ beak ged rc aal ‘tt 
a ae, ; it 


MOL t37915 OST o' se ay ia Ea Co i = rig 


On: by Te . Lays) oe Boal 
= 


a 


at [oval saol sri rg BA seed ginbafoo. to 4 
i 
tales ot bas (Sv! Stusges of Uae £5955 et at 23tde lalts vi af 


SUL Y iJ oo ; rag, Fons tol Boe sc acai ash 30 


omdebh soneckhss avec 1 2, a la Loeved set gakatase 3 ef a 183 


sian Sr wom Bannecs sits ont reakbiedar dase italia aed pie’ os 
> a oe 
tu500 mS ugil vd. bédayleps it accinecnaiaule if 
: ai, 
oom? us Tovensriw 990 Liiw sesll .fsval wt eSetr nohbee vhs 
' a 
m 


.beresnid@ocs SEB atisev Yitet suetsve Svods 16 asitea SB ‘x0 Vasey th 
. 
vbsem s4 wolle? Isysl saat onragad od to zautev 


2 bas ry 
fiiant aaptepood wet 2 no stardsh ah tud etdsaogso sonup cota beats: 


i cas Wipe, 
ited bidew eveavire bac abvours: sTok afsved % aa $5 od% wos} * 
7 - th ts i 
—o : oe wry 


~pivat sdf «.isbom oft Im t#63 2Si3od « ovite be 


67 


0861 


OL6T 


O96) 


(S9°0 = d) STeASeT ayeT 19q0790 aye] BUTY00) :0T TNNDIA 


ULYO IWIIYOLSIH ONISN S13A371 3yyT 


YbsA 


O€6T 0261 OT6I 


OS6T Ovét 068T 088 


O06 


+----+ Todd] Pepz0odv9y 


ToAsT pe ,eUTIAsyY 


9EL vel EL O&L 
T3A31 SMU 


SEL 


Ord 


oy ¢ 6 [oat : 


» en i 


aso. 1830 1a¢0 raep 


F870 


fEKE TEAET2 eine ut2ioKlceur Dvle 
AEG 


ErGHvE IQ? cocks srs Oem eee elle 


1300 


a0 © Mate 


| 
a ad a4 & ite Gab. Oe aw 
¢& | mo 
3 ‘7% i? 9 } aa vee ad = : 


oa F ve ; ye Pu = _ nF 


- sad ieee - 7" eat - aoe : 7 
a val o ne ‘ tel a vi? j 


—< 
— 
i?) 
iy 


>< 


an 


OL61 


0967 


OS6I 


(4 O[QeLLTALA) A9qGORDQ AOF S]oADY oyeY Buryoog 


YbsA 


> EL 


Ov6t O€6I 0c6t OT6t 


HUnoLel 


OO6T 


0681 


0881 


BEL g 
Sut Shel 


Ove 


CL 


qu 


ii ofdnbxyev) cvodeash wt eloved,clnd getsoad :44 ARTS — 


69 


tions between estimated and measured lake levels can be partially 
explained if the physical characteristics of the basin are investigated. 
Such an investigation will undoubtedly reveal numerous cause and effect 
relationships and at least several hypotheses could be postulated which 
could adequately explain the difference in model predicted lake levels 
and those recorded. Basic to most of the hypotheses that could be 
postulated would be the following observations. 

(1) The nature of the model is such that errors in the 
evapotranspiration estimation (hence runoff) for a given 
time period are ignored. This is due to the unavailability 
of the essential hydrological and meteorological information 
necessary to produce a model which can estimate runoff 
and model the physical properties of the basin correctly. 
The assumption of an average capacity of 10 cms may not 
be valid for most of the sub-basins in the Moraine. It 
is also generally accepted that Thornthwaite's procedure 
gives an over estimate of evapotranspiration therefore under 
estimating runoff. It would therefore be expected that 
higher lake levels could be predicted if a better method 
of estimating evapotranspiration becomes available. 

(2) The effective drainage area factor, F, cannot be thought 
of as a constant, for it obviously varies from one year to 
the next depending upon the relative amounts of precipitation 
and evapotranspiration that have occurred. This factor 
would also be expected to have a yearly variation with a 
maximum in the spring and a minimum in the winter when 


no runoff occurs at all. Within the summer rainfall period 
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this factor is also believed to have a significant 
variation. This summer variation being dictated by the 
duration, intensity and spatial arrangement of precipi- 
tation within the basin. 

(3) The dependence of runoff on the amount and more signifi- 
cantly, with respect to lake level changes, the location 
of rainfall within the watershed will also have its effect 
on the magnitude of lake level fluctuations. The spatial 
variation of rainfall and the topographic features of the 
basin such as areas of bogs, local depressions, drainage 
channels and outcrops will combine to further help in 
either the over or under-estimation of runoff amounts. 

(4) The amount of snowfall within the basin is a random variable 
which has been modeled by showing, with a suitable 
transformation, that it follows the normal distribution. 
Snow has however the unfortunate property of drifting into 
coulees and depressions whenever a wind blows. These 
deep banks localize runoff so that they yield concentrated 
pockets of snowmelt runoff. The extent of drifting will 
obviously vary from year to year depending upon the wind 
pattern and time of snowfall. This introduces another 
degree of unpredictability into the model and can only be 
assimulated into the water-balance model by physical 
measurement either directly by snow surveys or indirectly 


by one of the photogrammetric techniques. 
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(5) Contributions, from other lakes, to the water-balance of 
a lake are controlled by estimates of maximum lake levels. 
These estimates were made from topographical maps with a 
large contour interval. Consequently the contributions 
from one lake to another will not be correctly represented. 
It is believed that the method used underestimates the 
contribution from overflow and so, particularly for primary 
lakes lower down in the system, higher lake levels should 
be expected. 

This discussion could obviously be extended to include many of the 
other hydrological and meteorological variables affecting the hydrology 
of the Cooking Lake Moraine. Such further elaboration on these topics 
is not thought to be necessary as it is accepted that the water-balance 
model developed here has limited suitability for accurate lake level 
estimates due to the limited data available. However, the fluctuations 
in lakes levels, in particular the changes due to the randomness of 
natural meteorological variables such as temperature and precipitation, 
can be illustrated using the model and data generated via the Monte Carlo 


Method. 
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Simulated Lake Levels 

Many physical phenomena have a random or stochastic variability 
with a deterministic (often dominant) component which can be exactly 
defined. In nature many of the deterministic components of physical 
phenomena, such as lake levels, are influenced by man's everday life 
and the style of such a life. In the Cooking Lake Moraine man has 
altered the physical properties of the watershed by clearing land, 
developing agricultural lands and recreation areas, and diverting water 
to places outside of the natural boundaries. Developments of this 
kind can obviously have compensating effects when viewed comparatively 
but they do undoubtedly change the deterministic component of the many 
variables effecting the hydrology of the area. 

The alteration of the measurable component of a random variable can 
be subtle, such as increasing use of fertilisers could increase nutrient 
levels in surface and sub-surface runoff, and so enhance eutrophication 
consequently changing the evaporation regime. Or, the alteration can 
be drastic and possibly irreversible such as gross water diversions or 
drainage of marshy lands. Such a drastic change in the deterministic 
component of the lake level of Miquelon Lake occurred in 1927 when 
Camrose built a canal, partially emptying the lake. 

It is not possible to know what future flows in rivers or what 
precipitation events will be, but it is probable that future events will 
have the same stochastic properties as the observed historical record. 
This assumption forms the basis of stochastic hydrology and it is from 


this premise that lake level fluctuations are simulated. This model is 
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not intended to predict lake levels of the future because the effects 

and occurrence of the many variables are random and so are unpredictable, 
but the model shows that large variations in lake levels can be explained 
from meteorological variability. 

Program, MASTER, is adaptable to producing a long record of potential 
evaporation, surplus and deficit amounts, and lake level estimates by 
including a call to GRAND (page 55). When the values of the mean and 
standard deviation of the transformed variable, total monthly precipitation 
and the mean and standard deviation of mean monthly temperature are fed 
into the model a series of monthly lake levels for a specified period 
(500 years) is produced. Due to a limitation on the plotting subroutine 
only 500 points can be plotted. This limits the plotted output to 
either, approximately forty years of month to month fluctuations or to 
a plot of the lake level in a particular month for 500 years. The latter 
format has been adopted. 

Using the Monte Carlo Method described in Chapter 4 a 500 year 
series of total monthly precipitation and mean monthly temperature was 
obtained. This generated series is then used as the input data for the 
water-balance model so that, estimates of monthly lake levels can be 
plotted for a given month. An inspection of the simulated lake levels 
obtained from this input of randomly generated meteorological data can 
now be made. The simulated lake levels for October are shown in Figure 
12 for Cooking Lake and Appendix 4 for the other major lakes. 

The simulated lake levels correspond with the historically estimated 
lake levels in several general features. The first similarity is the 


Occurrence of sudden increases in lake levels spanning a relatively short 


period of several years. These sudden fluctuations are similar to the 
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observed rise in the lake level recorded at the turn of this century. 
However, these rises are not of the same magnitude nor of the same 
magnitude as the observed rise. Over the five hundred year period 
Simulated, three of these sudden restorations in lake levels occurred. 

The first occurred about year number 165 (Figure 12 ) and this rise in 
level culminated a relatively stable lake level period which had last 
approximately fifty years. The second and third rises however, occurred 
after the lake levels had been in a declining phase for approximately 
forty years in the first case and eighty years in the second. The second 
and third sudden rises occurred approximately at year 280 and 410 
respectively. Positive fluctuations in level did occur during the de- 
clining phases mentioned above but they did not effectively cause a 
reversal in the declining trend. It would appear from the inferences made 
above that a declining trend can only be reversed or halted by the occurrence 
of an exceptional hydrological event such as a very wet year or a series 
of wet years. 

After each restoration of the lake level by heavy runoff the system 
seems to decline towards an equilibrium level. This can be observed for 
both the historical plot of lake levels figure 10 and the simulated plot, 
figure 12 . When the lake system reaches its apparent equilibrium level 
it fluctuates randomly for a varying time period until another exceptional 
rainy period occurs. The system does not appear to have a definite trend 
towards either a high lake level or to a low lake level. Large fluctuations, 
both positive and negative, can be explained by the occurrence of 
exceptional hydrological events such as droughts and floods. 


The decline of the lake system contained within the boundaries of 


bitonae we sbeoose git int etney cite PH RE = on | 
Ole bue UGS aeey ts yletemixoryde ss17r990 agen robbie t : Ran 

Sb aid garb wisso bbb fayal at. acre ous aeot 

8 ah yloviiya2?s ton bib wens we see: Sa re 4 

a Bent encuhesesienee. sid) mori teogqe bivow a0 wheat ghbabissb ae +2 
sonersnope sais a bated $0 baatave7 sd tmp ds basxd asantinh = 
aside 670 Tasy Jaw yay s be dope snare! a Lento sesor® 


sot bevetedo od mo ehtT .tavel mobrditiope ae po satioeb 02 | ) 
Jolla basalumte sf2 bre OL sind alors. del to aolq fantroseid orld 
lawent mpbrdsliupe dgetaiqs 272i Sorssey, nsveve sist sta gedW . at . 
en ea asdtouse {hiny bobisq Sid goriaey & 203, Uhnohnes . 
has<d tkahiab & sved od sioqqe Jon ascb. madeua SMD \jaaweno baka 
2nok rand onl Deeg -fevel set wel ¢ oa 30 Loum et igh 8 redake 3 


Ye spmatawooe agld yt. eka BD Dome 
2booii bos -estlasoib ea 


76 


the Cooking Lake Moraine, over the last seventy years, can be explained, 
not from an investigation of land-use changes but from a more basic 
explanation based on the variability of natural phenomena. The simulated 
lake levels provide some evidence for a belief that the lake system 
prefers to decline from a high lake level to an equilibrium level around 
which it will oscillate until interrupted by an extreme wet year or wet 
years. This extreme hydrological event restores the lake system to 

an unstable high level from which the cycle starts again. However, this 
cycle cannot be assigned a certain periodicity as the culmination of a 
declining stage or equilibrium stage is based solely on the occurrence 


of a random event. 
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CHAPTER 6 


CONCLUSIONS AND RECOMMENDATIONS 


Conclusions 

The water-balance model developed for the Cooking Lake Moraine, shows 
the response of a lake system (with assumed steady-state physical 
properties) when excited by mathematically generated random meteorological 
Podabliest The response of such a system to a random input has been shown 
to be governed primarily by the variability of the natural random variables 
of precipitation and temperature. As a consequence of this dependence 
of the lake levels on the natural variability of climate, the response 
of the system is also random. 

An investigation was made to determine the frequency distribution of 
the random variables, total monthly precipitation and mean monthly temper- 
ature. It was shown that after determined transformations, the distribution 
of the two random variables can be fitted to the normal distribution. 

The characteristics of the simulated lake levels are similar in 
general features to a series of lake levels obtained by using historical 
data as input to a water-balance model. Lake levels appear to follow a 
cycle starting with a high unstable level. This unstable state is followed 
by a declining phase which eventually leads to an oscillatory stage where 
increases in level roughly balance the decreases. The latter two stages 
are subject to interruption, at any time, by an exceptional hydrological 
event. This event restores the lake level to a high unstable state and 
so effectively restarting the cycle. However no definable period can 
be assigned to this cycle as the exceptional hydrological event can occur 


at any time. 
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The effects of land use changes are not thought to play an important 
role in lake level fluctuations. However, land use changes affect the 
physical properties of the watershed and so the deterministic component 
of the water-balance model. Land use changes may affect the long term 
condition of the lake system by altering the evapotranspiration regime of 
the watershed. The water-balance components of precipitation, evapo- 
transpiration and runoff have the more significant effect on lake level 
fluctuations. The natural variability of precipitation and temperature 
impose on the lake level fluctuations a character which is more changeable 


or random than would be expected solely from land use alterations. 


Recommendations 

It is necessary in any future study to develop more refined methods 
in measuring the physical characteristics of the watershed. The properties 
of the basin which should be better defined include, (a) the areal extent 
of differing vegetal groups and soil types, (b) areas of depressional 
storage which do not contribute to lake runoff, and (c) the accurate 
survey of outlet control levels of each lake. If this elaboration was 
achieved, the water-balance model would be greatly improved with respect 
to both accuracy and reliability. 

A better evaporation estimation procedure is available in the Penman 
Method but it is believed that if a better model of the physical character 
of the watershed could be developed, the Thornthwaite procedure would still 
yield valuable results. 

The Monte Carlo Method has been used to generate normally distributed 
psuedo-random variables with a given mean and standard deviation. As it 


is possible that significant correlation could exist between meteorological 
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variables measured in consecutive months it would be advantageous to 
investigate this possibility in any future studies. If some correlation 
did exist to some considerable level the Monte Carlo Method should contain 


the lag-1 serial correlation coefficient. 
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Main Computer Program — MASTER 


DIMENSION PPTBAR(12) ,STPET(12) ,PPT(12,5€0) , TEMP(12,500) , 
*AHI (500) ,TEST (12) ,HRS(12) ,HI (12) ,A (500) ,E (12,500) ,PE (12,500), 
*SURP (12,500) ,SC (12,500) , DEF (12,569) ,STOR(13,500) , BET (12) , 
*SCHLL (500) ,CHLL(12,5¢0) ,AL(12,5C0) ,B(500) ,X (501) , ALPHA(29), 
*ALMAX (6) ,ALMIN (6) ,C(6),REST(6),TECA (6) ,ALOW(6),TEBAR (12), 
*AUP (6) ,1(509),P(5C9) , BLOW (6) , BHI (6) 

READ NUMEER CF YEARS GENERATED, AND THE VALUES OF THE MEAN 
AND STANDARD LEVIATION CF TOTAL MONTHLY TRANSFORMED PPT, 
REAL (5,499) N,MON 
499 FCRMAT (215) 
REAL (5,505) (ALPHA (I) ,1=1, 20) 
505 FCRMAT(20A4) 
READ (5,5(0) (PPTBAR (J) ,J=1,12) 
READ (5,560) (STPPT (J) ,J=1,12) 
500 FCRMAT(1ZF8.4) 
GENERATE N YRS CF NORMAILY DISTRIBUTED VALUES OF TRANSFORMED PPT, 

IX=6£539 

DO 1€9 I=1,12 

DC 101 J=1,N 

CALL GAUSS (IX,STPPT (I) ,PETBAR(I) ,P(J)) 

101 PET (1,J)=P(J) 

100 CONTINUE 

READ MEAN AND STANDARD DEVIATICN OF MEAN MONTHLY TEMPERATURE 
READ(5,50C) (TEBAR(I) ,1=1,12) 
READ (5,500) (TEST (I) ,1=1, 12) 

GENERATE N YRS OF MEAN MONTHLY TEMPERATURE 

1X=524287 

DEP 192 OT 4172 

DC 103 J=1,N 

CALL GAUSS (1IX,TEST (I) ,TEEAR(I) ,T(J)) 

103 TEME (I,J) =T (J) 

102 CCNIINUE 

READ IN KUMBER CF POSSIBLE HCURS OF SUNSHINE FOR EACH MONTH AT 5¢ 
READ (5,504) HRS 

504 FORMAT (12F6. 2) 
pe 104 J=1,N 
AHI (J) =€.0 

DETERPINE ANNUAL HEAT INDEX 

D619 5°12=1/)12 

TEMP (I,J) = (TEMP (I,J) -320) *50/9e 

IF (TEMPE (I,J) »LE.9.9) GOTC10 

HI (1) = (TEMP (I,J) /5e) **1651 

GCTO105 

190 HI(I)=0.0 
1405 AHI (J) =AHI(J) +HI(1) 
CTO=0.0179¥AHI (J) +9.492 
A (J) = (67 SE-08) *AHI (J) **3- (77, 1E-06) * AHI (J) **2+CTO 
ESTIMATE MONTHLY POTENTIAL EVAEORATION BY THORNTHWAITE'S METHCD 

p60106 2724312 

IF (TEMP (I,J) »LE.0.0) GOTO11 

B(I,d)=1.6*(10.9*TEME (I,J) /AHI (J) ) **A (J) 

PE (1,0) =HRS (I) *£(I,J) 

GC1C106 

44. PHYT)3)=0.0 
196 CCNTINUE 
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104 CCNTINUE 94 

WRITE MONTHLY ESTIMATES CF PE 
WRITE (6,6€9) 

600 FCRMAT(1X,*PCTENTENTIAL EVAPCTRANSPIRATICN*®/3X) 
WRITE (6,601) EE 

601 FCRMAT(1X,12F8. 2) 

REAL (5,500) BET 

READ IN MONTHLY TOTAL PET. THANSFORMATICNS 

TRANSFCRM GENERATED PPT. VALUES TO NATURAL VALUES 
BG IGG? OJ=1,N 
DG 108 I=1,12 
IF (FET (I,J) .18.0.0) GCTO12 
BET (1) =1.0/BET (1) 

PPT (1,3) =FPT (I,J) **BET (1) *2.54 
GCTC108 
12 PPI (I,J)=0.9 

108 CCNTINUE 

107 CCNTINUE 

ESTIMATE CHANGE IN STORAGE 
DC 149 J=1,N 
DG 110 IT=1,12 

119 SC (I,J) =PPT(1,J) -PE(I,3) 

109 CCNTIINUE 

SET STORAGE CAPACITY CF GHOTINED AND STORAGE CONDITION AT START 
CCVER=16.0 
AST=19,9 

ESTIMAT SURELUS OR DEFICIT FOR EACH MCNTH 
DC 1117 J=1,N 
S108 (1,J) =COVER 
DOU) PAA Pe 
STITCH {I+1,J) =STOR (I,J) +SC (I,J) 
LESIOM E+, 0) 936 777 

5 DEF (1,3) =-1.0*STOR (I+1,J) 
Sua (10) =0.0 
STOR ({I+1,0)=0.9 
GOIO112 
6 SURE (I,J)=9.9 
DEF (1,J)=9.9 
GCIC112 
7 IF(STOR(1+1,J) .LT.AST) GOTO6 
SUse4 D0) =STOR (Ue J) -AST 
DEF (I,J)=0.9 
STOR (1+1,J) =AST 
112 CCNTIINUE 
111 COVER=SiTCR (13,4) 
WRITE (6,498) N 

498 FCRMAT(1X,15) 

WRITE ESTIMATES OF SURPLUS AND DEFICIT 
WRITE (6,501) ((SC(I,J) ,1=1,12) ,J=1,N) 
WRITE (6,501) ((SURP (I,J) ,1=1,12) ,J=1,N) 

501 FORMAT (1X, 1258. 3) 

X(1)=9.9 
pO 113 I=1,N 
413 Meee 1) =X (TT) +1.0 


DETERMINE CHANGES IN LAKE LEVEL FOR GENERATED VALUES 
DC 120 K=1,6 
D=3.9 
NF=1 


READ (5,516) TELA (K) 
516 FCRMAT(F10.¥) 
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116 REAL (5,515) ALMAX {K) , ALMIN (K) ,C (K) , REST (K) a2 


Agi REACT (5,515) ALOW (K) ,AUP(K) , BLOW (K) , BHI (K) 
1718 515 FCRMAT (4F15. 4) 

119 DC 114 J=1,N 

120 DO 119) I=4,42 

121 AL (1,3) =C {K) +SuURP (I,J) *D/100.0 +SC (I,J) /100.0 
122 IF (AL(I,d) »GE.~ALMAX (K) ) GCTC239 
Zo IF (AL (I,J) LEeALMIN (K) )GCTC204 
124 CHLI (I,J) =AL (I,J) -C (K) 

125 GCici15 

126 230 CHLI (I,J) =ALMAX (K) -C (K) 

127 AL (I,J) =ALMAX (Kk) 

128 GCTC115 

129 204 CHILI (I,J) =).4 

130 AL (1,J) =ALMIN (K) 

131 115 C(K)=AL (I,J) 

2 Sar IF {C (K) »GE. REST (K) ) GOTO232 

133 SA=AICh (K) *C (K) ~BLOW (K) 

134 GCTO323 

135 232 SA=AUEF {K) *C (K) -BHIT (K) 

136 323 D=(TEDA(K)-SA) /SA 

137 4119 CCNIINUE 

730 114 CCNTINUE 

439 DC 116 J=1,N 

140 SCHLL (J) =0.0 

141 DEM?) £21, 42 

142 177 SCHII (J) =SCHLL (JpFCHLi (i, J) 

143 116 CCNTINUE 

144 DC 118 J=1,N 

145 118 B({J)=AL (MCN,J) 

146 C PLOT LAKE LEVELS FOR A PARTICULAR MONTH OF THE YEAR 
4147 CALL -CGEL (X,B,X,N,NF,1,1,4,1,0.0,50.0,10.0, 
148 * ALMIN(K) ,2.0,5.0,ALEHA, 6) 
149 120 CCNIINUE 

156 CALL CGPL (X, By XpNoO9 1414 4517,00 5 D0 e904 90,00,9.,ALPHA, 6) 
s Rese STOCE 

1&2 END 

153 SUBRCUTINE GAUSS (IX,5,AM,V) 

154 A=9.¢ 

455 De 59 I=1,12 

156 CALL KANDU(IX,1IY,Y) 

ibe Rf IK=I1Y 

158 506 A=A+tY 

159 V=(A-6.0) *St+AM 

160 RETUBN 

161 END 

162 SUCBRCUTINE RANDU(IX,IY,YFL) 

163 TY=IX*65539 

164 IF (IY)5,6,6 

165 5 LTY=LY+2147483647+1 

166 6 YFL=FLCAT (TY) 

167 YFL=Y¥FL*. 4656613E-9 

168 RETURN 

169 END 


END CF FILE 
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LAKE LEVELS OF OTHER LAKES IN THE SYSTEM 
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